library(raster)
## Loading required package: sp
library(tidyverse)
## ── Attaching packages ──────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.2     ✓ purrr   0.3.4
## ✓ tibble  3.0.3     ✓ dplyr   1.0.1
## ✓ tidyr   1.1.1     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0
## ── Conflicts ─────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x tidyr::extract() masks raster::extract()
## x dplyr::filter()  masks stats::filter()
## x dplyr::lag()     masks stats::lag()
## x dplyr::select()  masks raster::select()
library(getlandsat) 
library(sf) 
## Linking to GEOS 3.8.1, GDAL 3.1.1, PROJ 6.3.1
library(mapview)
## GDAL version >= 3.1.0 | setting mapviewOptions(fgb = TRUE)
library(osmdata)
## Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright
library(factoextra)
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
library(leaflet)

#Question 1

bb = read_csv('../data/uscities.csv') %>% 
  filter(city == "Palo") %>% 
  st_as_sf(coords = c("lng", "lat"), crs = 4326) %>% 
  st_transform(5070) %>% 
  st_buffer(5000) %>% 
  st_bbox() %>% 
  st_as_sfc() %>% 
  st_as_sf()
## Parsed with column specification:
## cols(
##   city = col_character(),
##   city_ascii = col_character(),
##   state_id = col_character(),
##   state_name = col_character(),
##   county_fips = col_double(),
##   county_name = col_character(),
##   county_fips_all = col_character(),
##   county_name_all = col_character(),
##   lat = col_double(),
##   lng = col_double(),
##   population = col_double(),
##   density = col_double(),
##   source = col_character(),
##   military = col_logical(),
##   incorporated = col_logical(),
##   timezone = col_character(),
##   ranking = col_double(),
##   zips = col_character(),
##   id = col_double()
## )

#Question 2

bbwgs = bb %>% st_transform(4326)
bb = st_bbox((bbwgs))

q2_image = getlandsat::lsat_scenes() %>% 
  filter(min_lat <= bb$ymin, min_lon <= bb$xmin) %>% 
  filter(max_lat >= bb$ymax, max_lon >= bb$xmax) %>% 
  filter(as.Date(acquisitionDate) == as.Date("2016-09-26"))
## Parsed with column specification:
## cols(
##   entityId = col_character(),
##   acquisitionDate = col_datetime(format = ""),
##   cloudCover = col_double(),
##   processingLevel = col_character(),
##   path = col_double(),
##   row = col_double(),
##   min_lat = col_double(),
##   min_lon = col_double(),
##   max_lat = col_double(),
##   max_lon = col_double(),
##   download_url = col_character()
## )
write.csv(q2_image, file = "../data/palo-flood-scene.csv")

Step 2

md = read_csv("../data/palo-flood-scene.csv")
## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   X1 = col_double(),
##   entityId = col_character(),
##   acquisitionDate = col_datetime(format = ""),
##   cloudCover = col_double(),
##   processingLevel = col_character(),
##   path = col_double(),
##   row = col_double(),
##   min_lat = col_double(),
##   min_lon = col_double(),
##   max_lat = col_double(),
##   max_lon = col_double(),
##   download_url = col_character()
## )
files = lsat_scene_files(md$download_url) %>% 
  filter(grepl(paste0("B", 1:6, ".TIF$", collapse= "|"), file)) %>% 
  arrange(file) %>% 
  pull(file)

st= sapply(files,lsat_image)
## File in cache
## File in cache
## File in cache
## File in cache
## File in cache
## File in cache
s=stack(st) %>% setNames(c(paste0("band", 1:6)))


cropper= bbwgs %>%  st_transform(crs(s))

r=crop(s, cropper)

#Question 3 True

par(mfrow=c(1,2))
plotRGB(r, r=4 ,g=3 ,b=2, stretch= "lin")
plotRGB(r, r=4 ,g=3 ,b=2, stretch= "hist")

Infrared

plotRGB(r, r=5 ,g=4 ,b=3, stretch= "lin")

plotRGB(r, r=5 ,g=4 ,b=3, stretch= "hist")

False Color

plotRGB(r, r=5 ,g=6 ,b=4, stretch= "lin")

plotRGB(r, r=5 ,g=6 ,b=4, stretch= "hist")

My Choice

plotRGB(r, r=2 ,g=4 ,b=6, stretch= "lin")

plotRGB(r, r=2 ,g=4 ,b=6, stretch= "hist")

Stretching determines the representation of each value. Histogram stretching tends to be lighter and brighter, while linear stretching seems to be darker

Question 4

5 Different Masks

palette= colorRampPalette(c("blue","white","red"))

ndvi=(r$band5 - r$band4) / (r$band5 + r$band4)
ndwi=(r$band3 - r$band5) / (r$band3 + r$band5)
mndwi=(r$band3 - r$band6) / (r$band3 + r$band6)
wri=(r$band3 + r$band4) / (r$band5 + r$band6)
swi= 1 / sqrt(r$band2 - r$band6)
## Warning in sqrt(getValues(x)): NaNs produced

Plotting

plot(ndvi,col=palette(256), legend=FALSE)

plot(ndwi,col=palette(256), legend=FALSE)

plot(mndwi,col=palette(256), legend=FALSE)

plot(wri,col=palette(256), legend=FALSE)

plot(swi,col=palette(256), legend=FALSE)

Water Thresholds

thresholding1= function(x){ifelse(x <= 0,1,NA)}
thresholding2= function(x){ifelse(x >= 0,1,NA)}
thresholding3= function(x){ifelse(x >= 0,1,NA)} 
thresholding4= function(x){ifelse(x >= 1,1,NA)}
thresholding5= function(x){ifelse(x <= 5,1,NA)}

flood1= calc(ndvi,thresholding1)
flood2= calc(ndwi,thresholding2)
flood3= calc(mndwi,thresholding3)
flood4= calc(wri,thresholding4)
flood5= calc(swi,thresholding5)

plot(flood1, col= "blue", legend=FALSE)

plot(flood2, col= "blue", legend=FALSE)

plot(flood3, col= "blue", legend=FALSE)

plot(flood4, col= "blue", legend=FALSE)

plot(flood5, col= "blue", legend=FALSE)

Question 5

Kmeans

set.seed(1)

r_extract = getValues(r)
r_nna = na.omit(r_extract) 


kmeans_r = kmeans(r_extract, 12)
## Warning: Quick-TRANSfer stage steps exceeded maximum (= 5882000)
fviz_cluster(kmeans_r, geom="point", data = r_extract)

thresholding_s= function(x){ifelse(x >= 0,1,0)}
flood_s= calc(ndwi,thresholding_s)
kmeans_raster = flood_s
values(kmeans_raster) = kmeans_r$cluster
plot(kmeans_raster, col = viridis::viridis(12))

Finding the correct KMeans cluster that represents our water pixels

com_table = table(values(flood_s), values(kmeans_raster))

which.max(com_table[2,])
## 12 
## 12
kmeans_raster[kmeans_raster != which.max(com_table[2,])] = 0
kmeans_raster[kmeans_raster != 0] =1
plot(kmeans_raster)

Changing remaining NAs to 0 for analysis

thresholdings1= function(x){ifelse(x <= 0,1,0)}
thresholdings3= function(x){ifelse(x >= 0,1,0)} 
thresholdings4= function(x){ifelse(x >= 1,1,0)}
thresholdings5= function(x){ifelse(x <= 5,1,0)}

floods1= calc(ndvi,thresholdings1)
floods3= calc(mndwi,thresholdings3)
floods4= calc(wri,thresholdings4)
floods5= calc(swi,thresholdings5)

fs = stack(floods1, flood_s, floods3, floods4, floods5, kmeans_raster)
plot(fs)

sum_fs = fs %>% sum() 
plot(sum_fs, col = blues9)

(cellStats(fs, sum) * res(fs)^2) /1e6
## layer.1 layer.2 layer.3 layer.4 layer.5 layer.6 
##  5.9994  6.4908 10.7451  7.6221 13.6809  7.9371

Extra Credit

How many layers classified the marked area in the video as a flood?

aoi = st_point(c(-91.78967, 42.06290)) %>% 
  st_sfc(crs = 4326) %>% 
  st_sf() %>% 
  st_transform(st_crs(fs))

raster::extract(sum_fs, aoi)
## [1] 3

There are 3 layers that classified this area as a flood; it is assumed that some rasters missed the marking of the area as a flood